library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(HDInterval) # for HPD intervals
library(posterior) # for posterior draws
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(ggeffects) # for partial effects plots
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # for the easystats ecosystem
library(ggridges) # for ridge plots
library(patchwork) # for multiple plots
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")Bayesian GLM Part6
1 Preparations
Load the necessary libraries
2 Scenario
An ecologist studying a rocky shore at Phillip Island, in southeastern Australia, was interested in how clumps of intertidal mussels are maintained (Quinn 1988). In particular, he wanted to know how densities of adult mussels affected recruitment of young individuals from the plankton. As with most marine invertebrates, recruitment is highly patchy in time, so he expected to find seasonal variation, and the interaction between season and density - whether effects of adult mussel density vary across seasons - was the aspect of most interest.
The data were collected from four seasons, and with two densities of adult mussels. The experiment consisted of clumps of adult mussels attached to the rocks. These clumps were then brought back to the laboratory, and the number of baby mussels recorded. There were 3-6 replicate clumps for each density and season combination.
| SEASON | DENSITY | RECRUITS | SQRTRECRUITS | GROUP |
|---|---|---|---|---|
| Spring | Low | 15 | 3.87 | SpringLow |
| .. | .. | .. | .. | .. |
| Spring | High | 11 | 3.32 | SpringHigh |
| .. | .. | .. | .. | .. |
| Summer | Low | 21 | 4.58 | SummerLow |
| .. | .. | .. | .. | .. |
| Summer | High | 34 | 5.83 | SummerHigh |
| .. | .. | .. | .. | .. |
| Autumn | Low | 14 | 3.74 | AutumnLow |
| .. | .. | .. | .. | .. |
| SEASON | Categorical listing of Season in which mussel clumps were collected independent variable |
| DENSITY | Categorical listing of the density of mussels within mussel clump independent variable |
| RECRUITS | The number of mussel recruits response variable |
| SQRTRECRUITS | Square root transformation of RECRUITS - needed to meet the test assumptions |
| GROUPS | Categorical listing of Season/Density combinations - used for checking ANOVA assumptions |
3 Read in the data
Rows: 42 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): SEASON, DENSITY, GROUP
dbl (2): RECRUITS, SQRTRECRUITS
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 42
Columns: 5
$ SEASON <chr> "Spring", "Spring", "Spring", "Spring", "Spring", "Spring…
$ DENSITY <chr> "Low", "Low", "Low", "Low", "Low", "High", "High", "High"…
$ RECRUITS <dbl> 15, 10, 13, 13, 5, 11, 10, 15, 10, 13, 1, 21, 31, 21, 18,…
$ SQRTRECRUITS <dbl> 3.872983, 3.162278, 3.605551, 3.605551, 2.236068, 3.31662…
$ GROUP <chr> "SpringLow", "SpringLow", "SpringLow", "SpringLow", "Spri…# A tibble: 6 × 5
SEASON DENSITY RECRUITS SQRTRECRUITS GROUP
<chr> <chr> <dbl> <dbl> <chr>
1 Spring Low 15 3.87 SpringLow
2 Spring Low 10 3.16 SpringLow
3 Spring Low 13 3.61 SpringLow
4 Spring Low 13 3.61 SpringLow
5 Spring Low 5 2.24 SpringLow
6 Spring High 11 3.32 SpringHighspc_tbl_ [42 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ SEASON : chr [1:42] "Spring" "Spring" "Spring" "Spring" ...
$ DENSITY : chr [1:42] "Low" "Low" "Low" "Low" ...
$ RECRUITS : num [1:42] 15 10 13 13 5 11 10 15 10 13 ...
$ SQRTRECRUITS: num [1:42] 3.87 3.16 3.61 3.61 2.24 ...
$ GROUP : chr [1:42] "SpringLow" "SpringLow" "SpringLow" "SpringLow" ...
- attr(*, "spec")=
.. cols(
.. SEASON = col_character(),
.. DENSITY = col_character(),
.. RECRUITS = col_double(),
.. SQRTRECRUITS = col_double(),
.. GROUP = col_character()
.. )
- attr(*, "problems")=<externalptr> quinn (42 rows and 5 variables, 5 shown)
ID | Name | Type | Missings | Values | N
---+--------------+-----------+----------+------------+-----------
1 | SEASON | character | 0 (0.0%) | Autumn | 10 (23.8%)
| | | | Spring | 11 (26.2%)
| | | | Summer | 12 (28.6%)
| | | | Winter | 9 (21.4%)
---+--------------+-----------+----------+------------+-----------
2 | DENSITY | character | 0 (0.0%) | High | 24 (57.1%)
| | | | Low | 18 (42.9%)
---+--------------+-----------+----------+------------+-----------
3 | RECRUITS | numeric | 0 (0.0%) | [0, 69] | 42
---+--------------+-----------+----------+------------+-----------
4 | SQRTRECRUITS | numeric | 0 (0.0%) | [0, 8.31] | 42
---+--------------+-----------+----------+------------+-----------
5 | GROUP | character | 0 (0.0%) | AutumnHigh | 6 (14.3%)
| | | | AutumnLow | 4 ( 9.5%)
| | | | SpringHigh | 6 (14.3%)
| | | | SpringLow | 5 (11.9%)
| | | | SummerHigh | 6 (14.3%)
| | | | SummerLow | 6 (14.3%)
| | | | WinterHigh | 6 (14.3%)
| | | | WinterLow | 3 ( 7.1%)
------------------------------------------------------------------| Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|
| RECRUITS | 25 | 0 | 18.3 | 15.7 | 0.0 | 13.5 | 69.0 |  |
| SQRTRECRUITS | 25 | 0 | 3.9 | 1.9 | 0.0 | 3.7 | 8.3 |  |
| N | % | |||||||
| SEASON | Autumn | 10 | 23.8 | |||||
| Spring | 11 | 26.2 | ||||||
| Summer | 12 | 28.6 | ||||||
| Winter | 9 | 21.4 | ||||||
| DENSITY | High | 24 | 57.1 | |||||
| Low | 18 | 42.9 | ||||||
| GROUP | AutumnHigh | 6 | 14.3 | |||||
| AutumnLow | 4 | 9.5 | ||||||
| SpringHigh | 6 | 14.3 | ||||||
| SpringLow | 5 | 11.9 | ||||||
| SummerHigh | 6 | 14.3 | ||||||
| SummerLow | 6 | 14.3 | ||||||
| WinterHigh | 6 | 14.3 | ||||||
| WinterLow | 3 | 7.1 |
quinn |>
dplyr::select(-SQRTRECRUITS) |>
modelsummary::datasummary_skim(
type = "numeric",
by = c("SEASON", "DENSITY")
)| SEASON | DENSITY | Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|---|---|
| RECRUITS | Spring | Low | 4 | 0 | 11.2 | 3.9 | 5.0 | 13.0 | 15.0 |  |
| Spring | High | 5 | 0 | 10.0 | 4.8 | 1.0 | 10.5 | 15.0 |  |
|
| Summer | Low | 5 | 0 | 22.0 | 6.1 | 14.0 | 21.0 | 31.0 |  |
|
| Summer | High | 6 | 0 | 48.2 | 15.0 | 28.0 | 51.5 | 69.0 |  |
|
| Autumn | Low | 4 | 0 | 18.2 | 3.1 | 14.0 | 19.0 | 21.0 |  |
|
| Autumn | High | 5 | 0 | 19.7 | 11.9 | 4.0 | 17.5 | 36.0 |  |
|
| Winter | Low | 2 | 0 | 2.7 | 4.6 | 0.0 | 0.0 | 8.0 |  |
|
| Winter | High | 5 | 0 | 5.7 | 3.3 | 1.0 | 5.0 | 10.0 |  |
4 Exploratory data analysis
Model formula: \[ \begin{align} y_i &\sim{} \mathcal{NB}(\lambda_i, \theta)\\ ln(\mu_i) &= \beta_0 + \sum_{j=1}^nT_{[i],j}.\beta_j\\ \beta_0 &\sim{} \mathcal{N}(2.4, 1.5)\\ \beta_{[1,2,3]} &\sim{} \mathcal{N}(0, 1)\\ \end{align} \]
where \(\beta_{0}\) is the y-intercept (mean of the first group), \(\beta_{[1,2,3]}\) are the vector of effects parameters (contrasting each group mean to that of the first group and \(T{[i],j}\) represents a \(i\) by \(j\) model matrix is a model matrix representing the season, density and their interaction on mussel recruitment.
5 Exploratory data analysis
The exploratory data analyses that we performed in the frequentist instalment of this example are equally valid here. That is, boxplots and/or violin plots for each population (substrate type).
# A tibble: 6 × 6
SEASON DENSITY RECRUITS SQRTRECRUITS GROUP fSEASON
<chr> <fct> <dbl> <dbl> <chr> <fct>
1 Spring Low 15 3.87 SpringLow Spring
2 Spring Low 10 3.16 SpringLow Spring
3 Spring Low 13 3.61 SpringLow Spring
4 Spring Low 13 3.61 SpringLow Spring
5 Spring Low 5 2.24 SpringLow Spring
6 Spring High 11 3.32 SpringHigh Spring 
Conclusions:
- there is clearly a relationship between mean and variance (as would be expected with the a Poisson
- evidently there are numerous zeros in the Winter/Low group # Fit the model
In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.
Priors for model 'quinn.rstanarmP'
------
Intercept (after predictors centered)
~ normal(location = 0, scale = 2.5)
Coefficients
Specified prior:
~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
Adjusted prior:
~ normal(location = [0,0,0,...], scale = [5.47,5.80,6.02,...])
------
See help('prior_summary.stanreg') for more detailsThis tells us:
for the intercept, when the family is Poisson, it is using a normal prior with a mean of 0 and a standard deviation of 2.5. The 2.5 is used for all intercepts. It is often scaled, but only if it is larger than 2.5 is the scaled version kept.
for the coefficients (in this case, the individual effects), the default prior is a normal prior centred around 0 with a standard deviations of 5.47, 5/8, 6.02 etc. This is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the numerical dummy variables for the predictor (then rounded).
fSEASONSummer fSEASONAutumn fSEASONWinter
5.467708 5.799380 6.019749
DENSITYLow fSEASONSummer:DENSITYLow fSEASONAutumn:DENSITYLow
4.991312 7.058781 8.414625
fSEASONWinter:DENSITYLow
9.590995 - there is no auxiliary prior as we are employing a Poisson distribution.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
quinn.rstanarm1 |>
ggpredict(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) |>
wrap_plots() &
scale_y_log10()Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## although, since there are zeros...
quinn.rstanarm1 |>
ggpredict(~ fSEASON + DENSITY) |>
plot(show_data = TRUE, jitter = FALSE) |>
wrap_plots() &
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
quinn.rstanarm1 |>
ggemmeans(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) |>
plot(show_data = TRUE) |>
wrap_plots() &
scale_y_log10()Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## although, since there are zeros...
quinn.rstanarm1 |>
ggemmeans(~ fSEASON + DENSITY) |>
plot(show_data = TRUE, jitter = FALSE) |>
wrap_plots() &
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Conclusions:
- we see that the range of predictions is fairly wide and the predicted means could range from 0 to very large (perhaps too large).
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
- \(\beta_0\): normal centred at 2.3 with a standard deviation of 5
- \(\beta_1\): normal centred at 0 with a standard deviation of 2
Remember the above are applied on the link scale.
I will also overlay the raw data for comparison.
quinn |>
group_by(fSEASON, DENSITY) |>
summarise(
Mean = log(mean(RECRUITS)),
SD = log(sd(RECRUITS))
)`summarise()` has grouped output by 'fSEASON'. You can override using the
`.groups` argument.# A tibble: 8 × 4
# Groups: fSEASON [4]
fSEASON DENSITY Mean SD
<fct> <fct> <dbl> <dbl>
1 Spring High 2.30 1.57
2 Spring Low 2.42 1.36
3 Summer High 3.87 2.71
4 Summer Low 3.09 1.81
5 Autumn High 2.98 2.48
6 Autumn Low 2.90 1.13
7 Winter High 1.73 1.20
8 Winter Low 0.981 1.53 (Intercept) fSEASONSummer fSEASONAutumn
Inf 6.024998 6.390475
fSEASONWinter DENSITYLow fSEASONSummer:DENSITYLow
6.633304 5.500045 7.778238
fSEASONAutumn:DENSITYLow fSEASONWinter:DENSITYLow
9.272276 10.568545 Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
quinn.rstanarm2 |>
ggpredict(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) |>
wrap_plots() &
scale_y_log10()Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## although, since there are zeros...
quinn.rstanarm2 |>
ggpredict(~ fSEASON + DENSITY) |>
plot(show_data = TRUE, jitter = FALSE) |>
wrap_plots() &
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
quinn.rstanarm2 |>
ggemmeans(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) |>
plot(show_data = TRUE) |>
wrap_plots() &
scale_y_log10()Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## although, since there are zeros...
quinn.rstanarm2 |>
ggemmeans(~ fSEASON + DENSITY) |>
plot(show_data = TRUE, jitter = FALSE) |>
wrap_plots() &
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Now lets refit, conditioning on the data.
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
prior class coef group resp dpar nlpar lb ub source
(flat) b default
(flat) b DENSITYLow (vectorized)
(flat) b fSEASONAutumn (vectorized)
(flat) b fSEASONAutumn:DENSITYLow (vectorized)
(flat) b fSEASONSummer (vectorized)
(flat) b fSEASONSummer:DENSITYLow (vectorized)
(flat) b fSEASONWinter (vectorized)
(flat) b fSEASONWinter:DENSITYLow (vectorized)
student_t(3, 2.6, 2.5) Intercept defaultRemember that the priors are applied on the link (in this case, log) scale.
quinn |>
group_by(fSEASON, DENSITY) |>
summarise(
Mean = mean(RECRUITS),
Median = median(RECRUITS),
MAD = mad(RECRUITS),
SD = sd(RECRUITS)
) |>
mutate(
log(Mean),
log(Median),
log(MAD),
log(SD)
)`summarise()` has grouped output by 'fSEASON'. You can override using the
`.groups` argument.# A tibble: 8 × 10
# Groups: fSEASON [4]
fSEASON DENSITY Mean Median MAD SD `log(Mean)` `log(Median)` `log(MAD)`
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Spring High 10 10.5 2.22 4.82 2.30 2.35 0.799
2 Spring Low 11.2 13 2.97 3.90 2.42 2.56 1.09
3 Summer High 48.2 51.5 15.6 15.0 3.87 3.94 2.75
4 Summer Low 22 21 6.67 6.13 3.09 3.04 1.90
5 Autumn High 19.7 17.5 12.6 11.9 2.98 2.86 2.53
6 Autumn Low 18.2 19 2.22 3.10 2.90 2.94 0.799
7 Winter High 5.67 5 3.71 3.33 1.73 1.61 1.31
8 Winter Low 2.67 0 0 4.62 0.981 -Inf -Inf
# ℹ 1 more variable: `log(SD)` <dbl>- \(\beta_0\): normal centred at 2.3 with a standard deviation of 1.5
- \(\beta_1\): normal centred at 0 with a standard deviation of 1
priors <- prior(normal(2.4, 1.5), class = "Intercept") +
prior(normal(0, 1), class = "b")
quinn.brm2 <- brm(quinn.form,
data = quinn,
prior = priors,
sample_prior = "only",
refresh = 0,
chains = 3,
iter = 5000,
thin = 5,
warmup = 2500,
backend = "rstan"
)Compiling Stan program...Start samplingData points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## Or since there are zeros
quinn.brm2 |>
ggpredict(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) +
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.Scale for y is already present.
Adding another scale for y, which will replace the existing scale.Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## Or since there are zeros
quinn.brm2 |>
ggemmeans(~ fSEASON + DENSITY) |>
plot(show_data = TRUE) +
scale_y_continuous(trans = scales::pseudo_log_trans())Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
[1] "b_Intercept" "b_fSEASONSummer"
[3] "b_fSEASONAutumn" "b_fSEASONWinter"
[5] "b_DENSITYLow" "b_fSEASONSummer:DENSITYLow"
[7] "b_fSEASONAutumn:DENSITYLow" "b_fSEASONWinter:DENSITYLow"
[9] "Intercept" "prior_Intercept"
[11] "prior_b" "lprior"
[13] "lp__" "accept_stat__"
[15] "stepsize__" "treedepth__"
[17] "n_leapfrog__" "divergent__"
[19] "energy__" 
quinn.brmP |>
posterior_samples() |>
dplyr::select(-`lp__`) |>
pivot_longer(everything(), names_to = "key") |>
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "b"),
Class = ifelse(str_detect(key, "Intercept"), "Intercept",
ifelse(str_detect(key, "b"), "b", "sigma")
),
Par = str_replace(key, "b_", "")
) |>
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge()) +
facet_wrap(~Class, scales = "free")Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
$N
[1] 42
$Y
[1] 15 10 13 13 5 11 10 15 10 13 1 21 31 21 18 14 27 34 49 69 55 28 54 14 18
[26] 20 21 4 22 30 36 13 13 8 0 0 10 1 5 9 4 5
$K
[1] 8
$Kc
[1] 7
$X
Intercept fSEASONSummer fSEASONAutumn fSEASONWinter DENSITYLow
1 1 0 0 0 1
2 1 0 0 0 1
3 1 0 0 0 1
4 1 0 0 0 1
5 1 0 0 0 1
6 1 0 0 0 0
7 1 0 0 0 0
8 1 0 0 0 0
9 1 0 0 0 0
10 1 0 0 0 0
11 1 0 0 0 0
12 1 1 0 0 1
13 1 1 0 0 1
14 1 1 0 0 1
15 1 1 0 0 1
16 1 1 0 0 1
17 1 1 0 0 1
18 1 1 0 0 0
19 1 1 0 0 0
20 1 1 0 0 0
21 1 1 0 0 0
22 1 1 0 0 0
23 1 1 0 0 0
24 1 0 1 0 1
25 1 0 1 0 1
26 1 0 1 0 1
27 1 0 1 0 1
28 1 0 1 0 0
29 1 0 1 0 0
30 1 0 1 0 0
31 1 0 1 0 0
32 1 0 1 0 0
33 1 0 1 0 0
34 1 0 0 1 1
35 1 0 0 1 1
36 1 0 0 1 1
37 1 0 0 1 0
38 1 0 0 1 0
39 1 0 0 1 0
40 1 0 0 1 0
41 1 0 0 1 0
42 1 0 0 1 0
fSEASONSummer:DENSITYLow fSEASONAutumn:DENSITYLow fSEASONWinter:DENSITYLow
1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
10 0 0 0
11 0 0 0
12 1 0 0
13 1 0 0
14 1 0 0
15 1 0 0
16 1 0 0
17 1 0 0
18 0 0 0
19 0 0 0
20 0 0 0
21 0 0 0
22 0 0 0
23 0 0 0
24 0 1 0
25 0 1 0
26 0 1 0
27 0 1 0
28 0 0 0
29 0 0 0
30 0 0 0
31 0 0 0
32 0 0 0
33 0 0 0
34 0 0 1
35 0 0 1
36 0 0 1
37 0 0 0
38 0 0 0
39 0 0 0
40 0 0 0
41 0 0 0
42 0 0 0
attr(,"assign")
[1] 0 1 1 1 2 3 3 3
attr(,"contrasts")
attr(,"contrasts")$fSEASON
Summer Autumn Winter
Spring 0 0 0
Summer 1 0 0
Autumn 0 1 0
Winter 0 0 1
attr(,"contrasts")$DENSITY
Low
High 0
Low 1
$prior_only
[1] 0
attr(,"class")
[1] "standata" "list" // generated with brms 2.22.0
functions {
}
data {
int<lower=1> N; // total number of observations
array[N] int Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> Kc; // number of population-level effects after centering
int prior_only; // should the likelihood be ignored?
}
transformed data {
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // regression coefficients
real Intercept; // temporary intercept for centered predictors
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += normal_lpdf(b | 0, 1);
lprior += normal_lpdf(Intercept | 2.4, 1.5);
}
model {
// likelihood including constants
if (!prior_only) {
target += poisson_log_glm_lpmf(Y | Xc, Intercept, b);
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// additionally sample draws from priors
real prior_b = normal_rng(0,1);
real prior_Intercept = normal_rng(2.4,1.5);
}6 MCMC sampling diagnostics
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
bayesplot MCMC module:
mcmc_acf
mcmc_acf_bar
mcmc_areas
mcmc_areas_ridges
mcmc_combo
mcmc_dens
mcmc_dens_chains
mcmc_dens_overlay
mcmc_hex
mcmc_hist
mcmc_hist_by_chain
mcmc_intervals
mcmc_neff
mcmc_neff_hist
mcmc_nuts_acceptance
mcmc_nuts_divergence
mcmc_nuts_energy
mcmc_nuts_stepsize
mcmc_nuts_treedepth
mcmc_pairs
mcmc_parcoord
mcmc_rank_ecdf
mcmc_rank_hist
mcmc_rank_overlay
mcmc_recover_hist
mcmc_recover_intervals
mcmc_recover_scatter
mcmc_rhat
mcmc_rhat_hist
mcmc_scatter
mcmc_trace
mcmc_trace_highlight
mcmc_violinOf these, we will focus on:
- mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Ratios all very high.
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Ratios all very high.
The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
always returns an ungrouped data frame and adjust accordingly.
ℹ The deprecated feature was likely used in the ggmcmc package.
Please report the issue at <https://github.com/xfim/ggmcmc/issues/>.
Ratios all very high.
The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
'pars' not specified. Showing first 10 parameters by default.
'pars' not specified. Showing first 10 parameters by default.
The chains appear well mixed and very similar
- stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
There is no evidence of auto-correlation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Ratios all very high.
7 Model validation
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
bayesplot PPC module:
ppc_bars
ppc_bars_grouped
ppc_boxplot
ppc_dens
ppc_dens_overlay
ppc_dens_overlay_grouped
ppc_ecdf_overlay
ppc_ecdf_overlay_grouped
ppc_error_binned
ppc_error_hist
ppc_error_hist_grouped
ppc_error_scatter
ppc_error_scatter_avg
ppc_error_scatter_avg_grouped
ppc_error_scatter_avg_vs_x
ppc_freqpoly
ppc_freqpoly_grouped
ppc_hist
ppc_intervals
ppc_intervals_grouped
ppc_km_overlay
ppc_km_overlay_grouped
ppc_loo_intervals
ppc_loo_pit_ecdf
ppc_loo_pit_overlay
ppc_loo_pit_qq
ppc_loo_ribbon
ppc_pit_ecdf
ppc_pit_ecdf_grouped
ppc_ribbon
ppc_ribbon_grouped
ppc_rootogram
ppc_scatter
ppc_scatter_avg
ppc_scatter_avg_grouped
ppc_stat
ppc_stat_2d
ppc_stat_freqpoly
ppc_stat_freqpoly_grouped
ppc_stat_grouped
ppc_violin_grouped- dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
The model draws are broadly similar to the observed data.
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.
The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments
- error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such.
- intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
The modelled predictions seem to do a reasonable job of representing the observations.
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
preds <- posterior_predict(quinn.rstanarm3, nsamples = 250, summary = FALSE)
quinn.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = quinn$RECRUITS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(quinn.resids)
Conclusions:
- the simulated residuals suggest there might be an issue of dispersion.
- it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
bayesplot PPC module:
ppc_bars
ppc_bars_grouped
ppc_boxplot
ppc_dens
ppc_dens_overlay
ppc_dens_overlay_grouped
ppc_ecdf_overlay
ppc_ecdf_overlay_grouped
ppc_error_binned
ppc_error_hist
ppc_error_hist_grouped
ppc_error_scatter
ppc_error_scatter_avg
ppc_error_scatter_avg_grouped
ppc_error_scatter_avg_vs_x
ppc_freqpoly
ppc_freqpoly_grouped
ppc_hist
ppc_intervals
ppc_intervals_grouped
ppc_km_overlay
ppc_km_overlay_grouped
ppc_loo_intervals
ppc_loo_pit_ecdf
ppc_loo_pit_overlay
ppc_loo_pit_qq
ppc_loo_ribbon
ppc_pit_ecdf
ppc_pit_ecdf_grouped
ppc_ribbon
ppc_ribbon_grouped
ppc_rootogram
ppc_scatter
ppc_scatter_avg
ppc_scatter_avg_grouped
ppc_stat
ppc_stat_2d
ppc_stat_freqpoly
ppc_stat_freqpoly_grouped
ppc_stat_grouped
ppc_violin_grouped- dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
The model draws appear to represent the shape of the observed data reasonably well
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.
Using all posterior draws for ppc type 'error_scatter_avg' by default.
The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments.
- intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
Using all posterior draws for ppc type 'intervals' by default.
Whilst the modelled predictions do a reasonable job of representing the observed data, the observed data do appear to be more varied than the model is representing.
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.quinn.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = quinn$RECRUITS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
quinn.resids |> plot()
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 2.7908, p-value < 2.2e-16
alternative hypothesis: two.sided
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 6.9444, p-value = 0.08
alternative hypothesis: two.sidedquinn.resids <- make_brms_dharma_res(quinn.brmP, integerResponse = TRUE)
wrap_elements(~ testUniformity(quinn.resids)) +
wrap_elements(~ plotResiduals(quinn.resids, form = factor(rep(1, nrow(quinn))))) +
wrap_elements(~ plotResiduals(quinn.resids, quantreg = TRUE)) +
wrap_elements(~ testDispersion(quinn.resids))
Conclusions:
- the simulated residuals do suggest that there might be a dispersion issue
- it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.
8 Explore negative binomial model
quinn.rstanarmNB <- stan_glm(RECRUITS ~ fSEASON * DENSITY,
data = quinn,
family = neg_binomial_2(link = "log"),
prior_intercept = normal(2.3, 2, autoscale = FALSE),
prior = normal(0, 10, autoscale = FALSE),
prior_aux = rstanarm::exponential(rate = 1, autoscale = FALSE),
prior_PD = FALSE,
refresh = 0,
chains = 3, iter = 5000, thin = 5, warmup = 2000
)
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
There seems to be a bug here. The expected values should be being back transformed to the response scale, however, they are clearly not. Notice that the expected values (and associated CI) are low and tiny respectively).
preds <- posterior_predict(quinn.rstanarmNB, nsamples = 250, summary = FALSE)
quinn.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = quinn$RECRUITS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(quinn.resids)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.25895, p-value = 0.03
alternative hypothesis: two.sidedNow possibly under-dispersed..
Warning: Found 2 observation(s) with a pareto_k > 0.7. We recommend calling 'loo' again with argument 'k_threshold = 0.7' in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 2 times to compute the ELPDs for the problematic observations directly.
Computed from 1800 by 42 log-likelihood matrix.
Estimate SE
elpd_loo -170.3 15.6
p_loo 23.6 5.4
looic 340.5 31.2
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.69] (good) 40 95.2% 143
(0.69, 1] (bad) 1 2.4% <NA>
(1, Inf) (very bad) 1 2.4% <NA>
See help('pareto-k-diagnostic') for details.Warning: Found 1 observation(s) with a pareto_k > 0.7. We recommend calling 'loo' again with argument 'k_threshold = 0.7' in order to calculate the ELPD without the assumption that these observations are negligible. This will refit the model 1 times to compute the ELPDs for the problematic observations directly.
Computed from 1800 by 42 log-likelihood matrix.
Estimate SE
elpd_loo -150.2 5.8
p_loo 8.9 3.2
looic 300.4 11.6
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.69] (good) 41 97.6% 672
(0.69, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 1 2.4% <NA>
See help('pareto-k-diagnostic') for details. elpd_diff se_diff
quinn.rstanarmNB 0.0 0.0
quinn.rstanarmP -20.1 12.0 quinn.form <- bf(RECRUITS ~ fSEASON * DENSITY, family = negbinomial(link = "log"))
get_prior(quinn.form, data = quinn) prior class coef group resp dpar
(flat) b
(flat) b DENSITYLow
(flat) b fSEASONAutumn
(flat) b fSEASONAutumn:DENSITYLow
(flat) b fSEASONSummer
(flat) b fSEASONSummer:DENSITYLow
(flat) b fSEASONWinter
(flat) b fSEASONWinter:DENSITYLow
student_t(3, 2.6, 2.5) Intercept
inv_gamma(0.4, 0.3) shape
nlpar lb ub source
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
0 defaultpriors <- prior(normal(2.4, 1.5), class = "Intercept") +
prior(normal(0, 1), class = "b") +
prior(gamma(0.01, 0.01), class = "shape")
quinn.brmsNB <- brm(quinn.form,
data = quinn,
prior = priors,
refresh = 0,
chains = 3,
iter = 5000,
thin = 5,
warmup = 2500,
backend = "rstan"
)Compiling Stan program...Start samplingWarning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.quinn.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = quinn$RECRUITS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(quinn.resids)
quinn.resids <- make_brms_dharma_res(quinn.brmsNB, integerResponse = TRUE)
wrap_elements(~ testUniformity(quinn.resids)) +
## wrap_elements(~plotResiduals(quinn.resids, form = factor(rep(1, nrow(quinn))))) +
wrap_elements(~ plotResiduals(quinn.resids, quantreg = TRUE)) +
wrap_elements(~ testDispersion(quinn.resids))
9 Partial effects plots
quinn.rstanarmNB |>
fitted_draws(newdata = quinn) |>
median_hdci() |>
ggplot(aes(x = fSEASON, colour = DENSITY, y = .value)) +
geom_pointrange(aes(ymin = .lower, ymax = .upper), position = position_dodge(width = 0.2)) +
geom_line(position = position_dodge(width = 0.2)) +
geom_point(data = quinn, aes(y = RECRUITS, x = fSEASON, colour = DENSITY), position = position_dodge(width = 0.2))Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
- Use [add_]epred_draws() to get the expectation of the posterior predictive.
- Use [add_]linpred_draws() to get the distribution of the linear predictor.
- For example, you used [add_]fitted_draws(..., scale = "response"), which
means you most likely want [add_]epred_draws(...).
NOTE: When updating to the new functions, note that the `model` parameter is now
named `object` and the `n` parameter is now named `ndraws`.
quinn.brmsNB |>
fitted_draws(newdata = quinn) |>
median_hdci() |>
ggplot(aes(x = fSEASON, colour = DENSITY, y = .value)) +
geom_pointrange(aes(ymin = .lower, ymax = .upper), position = position_dodge(width = 0.2)) +
geom_line(position = position_dodge(width = 0.2)) +
geom_point(data = quinn, aes(y = RECRUITS, x = fSEASON, colour = DENSITY), position = position_dodge(width = 0.2))Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
- Use [add_]epred_draws() to get the expectation of the posterior predictive.
- Use [add_]linpred_draws() to get the distribution of the linear predictor.
- For example, you used [add_]fitted_draws(..., scale = "response"), which
means you most likely want [add_]epred_draws(...).
NOTE: When updating to the new functions, note that the `model` parameter is now
named `object` and the `n` parameter is now named `ndraws`.
10 Model investigation
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
Model Info:
function: stan_glm
family: neg_binomial_2 [log]
formula: RECRUITS ~ fSEASON * DENSITY
algorithm: sampling
sample: 1800 (posterior sample size)
priors: see help('prior_summary')
observations: 42
predictors: 8
Estimates:
mean sd 10% 50% 90%
(Intercept) 2.3 0.3 2.0 2.3 2.6
fSEASONSummer 1.6 0.3 1.1 1.6 2.0
fSEASONAutumn 0.7 0.3 0.3 0.7 1.1
fSEASONWinter -0.6 0.4 -1.1 -0.6 -0.1
DENSITYLow 0.1 0.4 -0.3 0.1 0.6
fSEASONSummer:DENSITYLow -0.9 0.5 -1.5 -0.9 -0.3
fSEASONAutumn:DENSITYLow -0.2 0.5 -0.8 -0.2 0.5
fSEASONWinter:DENSITYLow -0.9 0.7 -1.7 -0.9 0.0
reciprocal_dispersion 3.9 1.2 2.5 3.8 5.6
Fit Diagnostics:
mean sd 10% 50% 90%
mean_PPD 19.3 3.1 15.6 18.9 23.3
The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
MCMC diagnostics
mcse Rhat n_eff
(Intercept) 0.0 1.0 1534
fSEASONSummer 0.0 1.0 1541
fSEASONAutumn 0.0 1.0 1657
fSEASONWinter 0.0 1.0 1560
DENSITYLow 0.0 1.0 1659
fSEASONSummer:DENSITYLow 0.0 1.0 1621
fSEASONAutumn:DENSITYLow 0.0 1.0 1741
fSEASONWinter:DENSITYLow 0.0 1.0 1654
reciprocal_dispersion 0.0 1.0 1652
mean_PPD 0.1 1.0 1538
log-posterior 0.1 1.0 1626
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Conclusions:
- the intercept represents the estimated mean of the first combination of Season (Spring) and Density (High). On the link scale this is 2.31
- the difference between Low and High adult density in spring is 1.59, although this is not significant
- the difference between Spring and Summer for High adult density is 0.69
- the difference between Spring and Autumn for High adult density is -0.58
- the difference between Spring and Winter for High adult density is 0.11
- if there were no interactions, we would expect the Low density Summer recruitment to be the additive of the main effects (Low and Summer). However, the modelled mean is 0.9 less than the additive effects would have expected. This value is significantly different to 0, indicating that there is evidence that the density effect in Summer is different to that in Spring.
- the density effect in Autumn and Winter were not found to be significantly different from what you would expect from an additive model.
tidyMCMC(quinn.rstanarmNB$stanfit, estimate.method = "median", conf.int = TRUE, conf.method = "HPDinterval", rhat = TRUE, ess = TRUE)# A tibble: 11 × 7
term estimate std.error conf.low conf.high rhat ess
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 (Intercept) 2.31 0.252 1.84 2.84 0.999 1534
2 fSEASONSummer 1.59 0.340 0.924 2.25 1.000 1541
3 fSEASONAutumn 0.686 0.342 0.0149 1.34 1.000 1657
4 fSEASONWinter -0.580 0.377 -1.30 0.203 0.998 1560
5 DENSITYLow 0.113 0.366 -0.595 0.833 0.999 1659
6 fSEASONSummer:DENSITYLow -0.899 0.496 -1.85 0.128 0.999 1621
7 fSEASONAutumn:DENSITYLow -0.176 0.522 -1.18 0.854 1.00 1741
8 fSEASONWinter:DENSITYLow -0.882 0.664 -2.04 0.648 0.999 1654
9 reciprocal_dispersion 3.93 1.22 1.73 6.17 1.00 1652
10 mean_PPD 19.3 3.15 13.8 26.0 1.00 1538
11 log-posterior -155. 2.41 -160. -151. 1.00 1626Conclusions:
See above
quinn.rstanarmNB$stanfit |>
summarise_draws(
median,
HDInterval::hdi,
rhat, length, ess_bulk, ess_tail
)# A tibble: 11 × 8
variable median lower upper rhat length ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.31 1.84e+0 2.84 0.999 1800 1526. 1408.
2 fSEASONSummer 1.59 9.24e-1 2.25 1.000 1800 1553. 1707.
3 fSEASONAutumn 0.686 1.49e-2 1.34 1.00 1800 1638. 1750.
4 fSEASONWinter -0.584 -1.30e+0 0.203 0.999 1800 1561. 1449.
5 DENSITYLow 0.107 -5.95e-1 0.833 0.999 1800 1623. 1665.
6 fSEASONSummer:DENS… -0.897 -1.85e+0 0.128 1.000 1800 1631. 1736.
7 fSEASONAutumn:DENS… -0.181 -1.18e+0 0.854 1.00 1800 1743. 1676.
8 fSEASONWinter:DENS… -0.895 -2.04e+0 0.648 1.00 1800 1636. 1542.
9 reciprocal_dispers… 3.76 1.73e+0 6.17 1.00 1800 1562. 1767.
10 mean_PPD 18.9 1.38e+1 26.0 1.00 1800 1572. 1784.
11 log-posterior -155. -1.60e+2 -151. 1.00 1800 1688. 1725.We can also alter the CI level.
quinn.rstanarmNB$stanfit |>
summarise_draws(
median,
~ HDInterval::hdi(.x, credMass = 0.9),
rhat, length, ess_bulk, ess_tail
)# A tibble: 11 × 8
variable median lower upper rhat length ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.31 1.88 2.71e+0 0.999 1800 1526. 1408.
2 fSEASONSummer 1.59 1.06 2.17e+0 1.000 1800 1553. 1707.
3 fSEASONAutumn 0.686 0.131 1.24e+0 1.00 1800 1638. 1750.
4 fSEASONWinter -0.584 -1.23 -1.12e-3 0.999 1800 1561. 1449.
5 DENSITYLow 0.107 -0.468 7.22e-1 0.999 1800 1623. 1665.
6 fSEASONSummer:DENS… -0.897 -1.77 -1.46e-1 1.000 1800 1631. 1736.
7 fSEASONAutumn:DENS… -0.181 -1.01 6.83e-1 1.00 1800 1743. 1676.
8 fSEASONWinter:DENS… -0.895 -2.02 1.59e-1 1.00 1800 1636. 1542.
9 reciprocal_dispers… 3.76 2.02 5.75e+0 1.00 1800 1562. 1767.
10 mean_PPD 18.9 14.4 2.42e+1 1.00 1800 1572. 1784.
11 log-posterior -155. -158. -1.51e+2 1.00 1800 1688. 1725.Arguably, it would be better to back-transform to the ratio scale
quinn.rstanarmNB$stanfit |>
summarise_draws(
~ median(exp(.x)),
~ HDInterval::hdi(exp(.x)),
rhat, length, ess_bulk, ess_tail
)# A tibble: 11 × 8
variable `~median(exp(.x))` lower upper rhat length ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Interce… 1.01e+ 1 5.82e+ 0 1.62e+ 1 0.999 1800 1526. 1408.
2 fSEASONS… 4.90e+ 0 2.22e+ 0 8.70e+ 0 1.000 1800 1553. 1707.
3 fSEASONA… 1.99e+ 0 8.83e- 1 3.58e+ 0 1.00 1800 1638. 1750.
4 fSEASONW… 5.58e- 1 1.92e- 1 1.04e+ 0 0.999 1800 1561. 1449.
5 DENSITYL… 1.11e+ 0 4.44e- 1 2.05e+ 0 0.999 1800 1623. 1665.
6 fSEASONS… 4.08e- 1 9.66e- 2 9.35e- 1 1.000 1800 1631. 1736.
7 fSEASONA… 8.34e- 1 1.96e- 1 1.99e+ 0 1.00 1800 1743. 1676.
8 fSEASONW… 4.09e- 1 5.38e- 2 1.30e+ 0 1.00 1800 1636. 1542.
9 reciproc… 4.30e+ 1 4.70e+ 0 4.38e+ 2 1.00 1800 1562. 1767.
10 mean_PPD 1.68e+ 8 4.69e+ 3 8.11e+10 1.00 1800 1572. 1784.
11 log-post… 6.06e-68 1.27e-75 1.01e-66 1.00 1800 1688. 1725.# A draws_df: 600 iterations, 3 chains, and 11 variables
(Intercept) fSEASONSummer fSEASONAutumn fSEASONWinter DENSITYLow
1 2.3 1.5 0.75 -0.43 -0.012
2 2.3 1.4 0.86 -0.32 0.439
3 1.9 2.1 1.09 -0.35 0.154
4 2.4 1.3 0.66 -0.80 -0.025
5 2.2 1.5 1.21 -0.70 -0.106
6 2.3 1.6 0.69 -0.45 0.062
7 2.4 1.6 0.53 -0.60 0.166
8 2.0 2.2 0.51 -0.61 0.868
9 2.2 1.9 0.78 -0.73 0.336
10 2.2 1.8 0.52 -0.33 0.269
fSEASONSummer:DENSITYLow fSEASONAutumn:DENSITYLow fSEASONWinter:DENSITYLow
1 -0.83 0.025 -0.34
2 -1.23 -0.651 -2.02
3 -0.80 -0.316 -1.34
4 -0.41 0.162 -0.34
5 -0.47 -0.480 -0.50
6 -0.75 0.073 -0.98
7 -1.55 -0.642 -0.81
8 -1.35 -0.721 -1.50
9 -1.79 -0.074 -1.02
10 -1.06 -0.300 -0.82
# ... with 1790 more draws, and 3 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}quinn.rstanarmNB$stanfit |>
as_draws_df() |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk
)# A tibble: 11 × 6
variable median lower upper rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.31 1.84 2.84 0.999 1526.
2 fSEASONSummer 1.59 0.924 2.25 1.000 1553.
3 fSEASONAutumn 0.686 0.0149 1.34 1.00 1638.
4 fSEASONWinter -0.584 -1.30 0.203 0.999 1561.
5 DENSITYLow 0.107 -0.595 0.833 0.999 1623.
6 fSEASONSummer:DENSITYLow -0.897 -1.85 0.128 1.000 1631.
7 fSEASONAutumn:DENSITYLow -0.181 -1.18 0.854 1.00 1743.
8 fSEASONWinter:DENSITYLow -0.895 -2.04 0.648 1.00 1636.
9 reciprocal_dispersion 3.76 1.73 6.17 1.00 1562.
10 mean_PPD 18.9 13.8 26.0 1.00 1572.
11 log-posterior -155. -160. -151. 1.00 1688.quinn.rstanarmNB$stanfit |>
as_draws_df() |>
exp() |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk
)# A tibble: 11 × 6
variable median lower upper rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.01e+ 1 5.82e+ 0 1.62e+ 1 0.999 1526.
2 fSEASONSummer 4.90e+ 0 2.22e+ 0 8.70e+ 0 1.000 1553.
3 fSEASONAutumn 1.99e+ 0 8.83e- 1 3.58e+ 0 1.00 1638.
4 fSEASONWinter 5.58e- 1 1.92e- 1 1.04e+ 0 0.999 1561.
5 DENSITYLow 1.11e+ 0 4.44e- 1 2.05e+ 0 0.999 1623.
6 fSEASONSummer:DENSITYLow 4.08e- 1 9.66e- 2 9.35e- 1 1.000 1631.
7 fSEASONAutumn:DENSITYLow 8.34e- 1 1.96e- 1 1.99e+ 0 1.00 1743.
8 fSEASONWinter:DENSITYLow 4.09e- 1 5.38e- 2 1.30e+ 0 1.000 1636.
9 reciprocal_dispersion 4.30e+ 1 4.70e+ 0 4.38e+ 2 1.00 1562.
10 mean_PPD 1.68e+ 8 4.69e+ 3 8.11e+10 1.00 1572.
11 log-posterior 6.06e-68 1.27e-75 1.01e-66 1.00 1688.Due to the presence of a log transform in the predictor, it is better to use the regex version.
[1] "(Intercept)" "fSEASONSummer"
[3] "fSEASONAutumn" "fSEASONWinter"
[5] "DENSITYLow" "fSEASONSummer:DENSITYLow"
[7] "fSEASONAutumn:DENSITYLow" "fSEASONWinter:DENSITYLow"
[9] "reciprocal_dispersion" "accept_stat__"
[11] "stepsize__" "treedepth__"
[13] "n_leapfrog__" "divergent__"
[15] "energy__" quinn.draw <- quinn.rstanarmNB |> gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE)
quinn.draw# A tibble: 14,400 × 5
# Groups: .variable [8]
.chain .iteration .draw .variable .value
<int> <int> <int> <chr> <dbl>
1 1 1 1 (Intercept) 2.33
2 1 2 2 (Intercept) 2.26
3 1 3 3 (Intercept) 1.94
4 1 4 4 (Intercept) 2.40
5 1 5 5 (Intercept) 2.24
6 1 6 6 (Intercept) 2.26
7 1 7 7 (Intercept) 2.42
8 1 8 8 (Intercept) 1.97
9 1 9 9 (Intercept) 2.19
10 1 10 10 (Intercept) 2.16
# ℹ 14,390 more rowsexceedP <- function(x, Val = 0) mean(x > Val)
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
summarise_draws(
median,
HDInterval::hdi,
rhat,
length,
ess_bulk,
ess_tail,
~ exceedP(.x, 1)
)# A tibble: 8 × 10
# Groups: .variable [8]
.variable variable median lower upper rhat length ess_bulk ess_tail
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) .value 10.1 5.82 16.2 0.999 1800 1526. 1408.
2 DENSITYLow .value 1.11 0.444 2.05 0.999 1800 1623. 1665.
3 fSEASONAutumn .value 1.99 0.883 3.58 1.00 1800 1638. 1750.
4 fSEASONAutumn:DE… .value 0.834 0.196 1.99 1.00 1800 1743. 1676.
5 fSEASONSummer .value 4.90 2.22 8.70 1.000 1800 1553. 1707.
6 fSEASONSummer:DE… .value 0.408 0.0966 0.935 1.000 1800 1631. 1736.
7 fSEASONWinter .value 0.558 0.192 1.04 0.999 1800 1561. 1449.
8 fSEASONWinter:DE… .value 0.409 0.0538 1.30 1.000 1800 1636. 1542.
# ℹ 1 more variable: `~exceedP(.x, 1)` <dbl>We can then summarise this
# A tibble: 8 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 2.31 1.84 2.84 0.95 median hdci
2 DENSITYLow 0.107 -0.595 0.833 0.95 median hdci
3 fSEASONAutumn 0.686 0.0187 1.35 0.95 median hdci
4 fSEASONAutumn:DENSITYLow -0.181 -1.18 0.854 0.95 median hdci
5 fSEASONSummer 1.59 0.924 2.25 0.95 median hdci
6 fSEASONSummer:DENSITYLow -0.897 -1.85 0.128 0.95 median hdci
7 fSEASONWinter -0.584 -1.28 0.222 0.95 median hdci
8 fSEASONWinter:DENSITYLow -0.895 -2.12 0.603 0.95 median hdci We could alternatively express the parameters on the response scale.
# A tibble: 8 × 7
.variable `exp(.value)` .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 10.1 5.69 16.0 0.95 median hdci
2 DENSITYLow 1.11 0.455 2.06 0.95 median hdci
3 fSEASONAutumn 1.99 0.967 3.73 0.95 median hdci
4 fSEASONAutumn:DENSITYLow 0.834 0.196 1.99 0.95 median hdci
5 fSEASONSummer 4.90 2.22 8.70 0.95 median hdci
6 fSEASONSummer:DENSITYLow 0.408 0.0966 0.935 0.95 median hdci
7 fSEASONWinter 0.558 0.192 1.04 0.95 median hdci
8 fSEASONWinter:DENSITYLow 0.409 0.0538 1.30 0.95 median hdci quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()
We could alternatively express the parameters on the response scale.
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
group_by(.variable) |>
mutate(.value = exp(.value)) |>
median_hdci()# A tibble: 8 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 10.1 5.69 16.0 0.95 median hdci
2 DENSITYLow 1.11 0.455 2.06 0.95 median hdci
3 fSEASONAutumn 1.99 0.967 3.73 0.95 median hdci
4 fSEASONAutumn:DENSITYLow 0.834 0.196 1.99 0.95 median hdci
5 fSEASONSummer 4.90 2.22 8.70 0.95 median hdci
6 fSEASONSummer:DENSITYLow 0.408 0.0966 0.935 0.95 median hdci
7 fSEASONWinter 0.558 0.192 1.04 0.95 median hdci
8 fSEASONWinter:DENSITYLow 0.409 0.0538 1.30 0.95 median hdci quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
ggplot() +
geom_vline(xintercept = 1, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
scale_x_continuous("", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
theme_classic()
## Link scale
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)Warning: `stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95), labels =
scales::percent_format()))` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()))` instead.
## Fractional scale
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
ggplot() +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
scale_x_continuous(trans = scales::log2_trans())
quinn.rstanarmNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
quinn.rstanarmNB |>
gather_draws(`.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
geom_vline(xintercept = 0, linetype = "dashed")
quinn.rstanarmNB |>
gather_draws(`.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = exp(.value), y = .variable)) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans())
quinn.rstanarmNB |>
gather_draws(`.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
geom_vline(xintercept = 0, linetype = "dashed")Picking joint bandwidth of 0.0876
## Or on a fractional scale
quinn.rstanarmNB |>
gather_draws(`.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_density_ridges_gradient(
aes(
x = exp(.value),
y = .variable,
fill = stat(x)
),
alpha = 0.4, colour = "white",
quantile_lines = TRUE,
quantiles = c(0.025, 0.975)
) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans()) +
scale_fill_viridis_c(option = "C")Picking joint bandwidth of 0.126
This is purely a graphical depiction on the posteriors.
# A tibble: 1,800 × 18
.chain .iteration .draw `(Intercept)` fSEASONSummer fSEASONAutumn
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 2.33 1.46 0.754
2 1 2 2 2.26 1.39 0.862
3 1 3 3 1.94 2.06 1.09
4 1 4 4 2.40 1.27 0.660
5 1 5 5 2.24 1.53 1.21
6 1 6 6 2.26 1.59 0.685
7 1 7 7 2.42 1.58 0.530
8 1 8 8 1.97 2.25 0.507
9 1 9 9 2.19 1.87 0.780
10 1 10 10 2.16 1.79 0.524
# ℹ 1,790 more rows
# ℹ 12 more variables: fSEASONWinter <dbl>, DENSITYLow <dbl>,
# `fSEASONSummer:DENSITYLow` <dbl>, `fSEASONAutumn:DENSITYLow` <dbl>,
# `fSEASONWinter:DENSITYLow` <dbl>, reciprocal_dispersion <dbl>,
# accept_stat__ <dbl>, stepsize__ <dbl>, treedepth__ <dbl>,
# n_leapfrog__ <dbl>, divergent__ <dbl>, energy__ <dbl># A tibble: 1,800 × 11
.chain .iteration .draw `(Intercept)` fSEASONSummer fSEASONAutumn
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 2.33 1.46 0.754
2 1 2 2 2.26 1.39 0.862
3 1 3 3 1.94 2.06 1.09
4 1 4 4 2.40 1.27 0.660
5 1 5 5 2.24 1.53 1.21
6 1 6 6 2.26 1.59 0.685
7 1 7 7 2.42 1.58 0.530
8 1 8 8 1.97 2.25 0.507
9 1 9 9 2.19 1.87 0.780
10 1 10 10 2.16 1.79 0.524
# ℹ 1,790 more rows
# ℹ 5 more variables: fSEASONWinter <dbl>, DENSITYLow <dbl>,
# `fSEASONSummer:DENSITYLow` <dbl>, `fSEASONAutumn:DENSITYLow` <dbl>,
# `fSEASONWinter:DENSITYLow` <dbl>Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.# A tibble: 1,800 × 9
`(Intercept)` fSEASONSummer fSEASONAutumn fSEASONWinter DENSITYLow
<dbl> <dbl> <dbl> <dbl> <dbl>
1 2.33 1.46 0.754 -0.427 -0.0120
2 2.26 1.39 0.862 -0.324 0.439
3 1.94 2.06 1.09 -0.349 0.154
4 2.40 1.27 0.660 -0.802 -0.0245
5 2.24 1.53 1.21 -0.703 -0.106
6 2.26 1.59 0.685 -0.450 0.0619
7 2.42 1.58 0.530 -0.601 0.166
8 1.97 2.25 0.507 -0.610 0.868
9 2.19 1.87 0.780 -0.729 0.336
10 2.16 1.79 0.524 -0.331 0.269
# ℹ 1,790 more rows
# ℹ 4 more variables: `fSEASONSummer:DENSITYLow` <dbl>,
# `fSEASONAutumn:DENSITYLow` <dbl>, `fSEASONWinter:DENSITYLow` <dbl>,
# reciprocal_dispersion <dbl>The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
Family: negbinomial
Links: mu = log; shape = identity
Formula: RECRUITS ~ fSEASON * DENSITY
Data: quinn (Number of observations: 42)
Draws: 3 chains, each with iter = 5000; warmup = 2500; thin = 5;
total post-warmup draws = 1500
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 2.42 0.19 2.06 2.82 1.00 1500
fSEASONSummer 1.41 0.25 0.91 1.89 1.00 1410
fSEASONAutumn 0.54 0.25 0.05 1.02 1.00 1487
fSEASONWinter -0.68 0.30 -1.28 -0.08 1.00 1510
DENSITYLow -0.06 0.26 -0.60 0.46 1.00 1502
fSEASONSummer:DENSITYLow -0.66 0.35 -1.33 0.05 1.00 1409
fSEASONAutumn:DENSITYLow 0.02 0.38 -0.69 0.75 1.00 1473
fSEASONWinter:DENSITYLow -0.61 0.48 -1.54 0.30 1.00 1553
Tail_ESS
Intercept 1347
fSEASONSummer 1461
fSEASONAutumn 1501
fSEASONWinter 1499
DENSITYLow 1440
fSEASONSummer:DENSITYLow 1439
fSEASONAutumn:DENSITYLow 1542
fSEASONWinter:DENSITYLow 1385
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 6.68 2.99 2.66 14.45 1.00 1422 1590
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
- the intercept represents the estimated mean of the first combination of Season (Spring) and Density (High). On the link scale this is 2.42
- the difference between Low and High adult density in spring is 1.41, although this is not significant
- the difference between Spring and Summer for High adult density is 0.54
- the difference between Spring and Autumn for High adult density is -0.68
- the difference between Spring and Winter for High adult density is -0.06
- if there were no interactions, we would expect the Low density Summer recruitment to be the additive of the main effects (Low and Summer). However, the modelled mean is 0.66 less than the additive effects would have expected. This value is significantly different to 0, indicating that there is evidence that the density effect in Summer is different to that in Spring.
- the density effect in Autumn and Winter were not found to be significantly different from what you would expect from an additive model.
quinn.brmsNB$fit |>
tidyMCMC(
estimate.method = "median",
conf.int = TRUE,
conf.method = "HPDinterval",
rhat = TRUE,
ess = TRUE
)# A tibble: 11 × 7
term estimate std.error conf.low conf.high rhat ess
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 b_Intercept 2.42 0.188 2.03 2.78 0.999 1475
2 b_fSEASONSummer 1.41 0.249 0.923 1.90 1.00 1392
3 b_fSEASONAutumn 0.541 0.247 0.0446 1.01 0.999 1482
4 b_fSEASONWinter -0.677 0.301 -1.22 -0.0439 0.999 1502
5 b_DENSITYLow -0.0584 0.262 -0.628 0.424 1.00 1502
6 b_fSEASONSummer:DENSITYLow -0.657 0.350 -1.34 0.0263 0.999 1321
7 b_fSEASONAutumn:DENSITYLow 0.0223 0.381 -0.699 0.744 1.00 1467
8 b_fSEASONWinter:DENSITYLow -0.609 0.482 -1.52 0.319 0.999 1550
9 shape 6.68 2.99 2.19 12.7 1.00 1479
10 Intercept 2.65 0.0803 2.49 2.80 0.999 1488
11 lprior -16.4 0.843 -18.1 -14.9 0.999 1376Conclusions:
see above
# A draws_df: 500 iterations, 3 chains, and 12 variables
b_Intercept b_fSEASONSummer b_fSEASONAutumn b_fSEASONWinter b_DENSITYLow
1 2.3 1.37 0.979 -0.36 -0.111
2 2.8 1.14 0.089 -0.96 -0.703
3 2.8 1.01 0.380 -0.56 -0.626
4 2.1 1.73 0.922 -0.47 0.414
5 2.2 1.57 0.697 -0.27 0.198
6 2.2 1.56 0.963 -0.68 0.266
7 2.5 1.17 0.572 -0.74 -0.209
8 2.8 1.26 0.625 -0.72 -0.324
9 2.8 0.78 0.754 -1.13 -0.596
10 2.5 1.28 0.509 -0.71 0.011
b_fSEASONSummer:DENSITYLow b_fSEASONAutumn:DENSITYLow
1 -0.252 -0.21
2 -0.172 0.59
3 0.109 0.35
4 -1.145 -0.40
5 -0.727 -0.25
6 -0.975 -0.21
7 -0.138 0.11
8 -0.690 0.31
9 0.041 0.11
10 -0.933 -0.26
b_fSEASONWinter:DENSITYLow
1 -0.172
2 -0.042
3 -0.825
4 -0.798
5 -1.097
6 -0.536
7 -0.106
8 -0.746
9 -0.610
10 -1.258
# ... with 1490 more draws, and 4 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}quinn.brmsNB |>
as_draws_df() |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk, ess_tail
)# A tibble: 12 × 7
variable median lower upper rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 2.42e+0 2.03e+0 2.78e+0 1.00 1500. 1347.
2 b_fSEASONSummer 1.42e+0 9.23e-1 1.90e+0 1.00 1410. 1461.
3 b_fSEASONAutumn 5.37e-1 4.46e-2 1.01e+0 1.00 1487. 1501.
4 b_fSEASONWinter -6.66e-1 -1.22e+0 -4.39e-2 1.00 1510. 1499.
5 b_DENSITYLow -5.13e-2 -6.28e-1 4.24e-1 1.00 1502. 1440.
6 b_fSEASONSummer:DENSITYLow -6.79e-1 -1.34e+0 2.63e-2 1.000 1408. 1439.
7 b_fSEASONAutumn:DENSITYLow 2.06e-2 -6.99e-1 7.44e-1 1.00 1473. 1542.
8 b_fSEASONWinter:DENSITYLow -6.06e-1 -1.52e+0 3.19e-1 0.999 1552. 1385.
9 shape 6.06e+0 2.19e+0 1.27e+1 1.00 1422. 1590.
10 Intercept 2.65e+0 2.49e+0 2.80e+0 1.00 1544. 1169.
11 lprior -1.64e+1 -1.81e+1 -1.49e+1 1.00 1398. 1576.
12 lp__ -1.57e+2 -1.62e+2 -1.53e+2 1.00 1475. 1500.quinn.brmsNB |>
as_draws_df() |>
exp() |>
summarise_draws(
median,
HDInterval::hdi,
rhat,
length,
ess_bulk, ess_tail
)# A tibble: 12 × 8
variable median lower upper rhat length ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 1.13e+ 1 7.63e+ 0 1.61e+ 1 1.00 1500 1500. 1347.
2 b_fSEASONSummer 4.13e+ 0 2.28e+ 0 6.29e+ 0 1.00 1500 1410. 1461.
3 b_fSEASONAutumn 1.71e+ 0 9.52e- 1 2.64e+ 0 1.00 1500 1487. 1501.
4 b_fSEASONWinter 5.14e- 1 2.31e- 1 8.34e- 1 1.00 1500 1510. 1499.
5 b_DENSITYLow 9.50e- 1 4.78e- 1 1.46e+ 0 1.00 1500 1502. 1440.
6 b_fSEASONSummer:DE… 5.07e- 1 2.29e- 1 9.57e- 1 1.000 1500 1408. 1439.
7 b_fSEASONAutumn:DE… 1.02e+ 0 4.46e- 1 2.02e+ 0 1.00 1500 1473. 1542.
8 b_fSEASONWinter:DE… 5.46e- 1 1.68e- 1 1.21e+ 0 0.999 1500 1552. 1385.
9 shape 4.27e+ 2 3.39e+ 0 2.27e+ 5 1.00 1500 1422. 1590.
10 Intercept 1.41e+ 1 1.21e+ 1 1.65e+ 1 1.00 1500 1544. 1169.
11 lprior 7.70e- 8 2.48e- 9 2.53e- 7 1.00 1500 1398. 1576.
12 lp__ 9.09e-69 8.07e-75 1.50e-67 1.00 1500 1475. 1500.Due to the presence of a log transform in the predictor, it is better to use the regex version.
[1] "b_Intercept" "b_fSEASONSummer"
[3] "b_fSEASONAutumn" "b_fSEASONWinter"
[5] "b_DENSITYLow" "b_fSEASONSummer:DENSITYLow"
[7] "b_fSEASONAutumn:DENSITYLow" "b_fSEASONWinter:DENSITYLow"
[9] "shape" "Intercept"
[11] "lprior" "lp__"
[13] "accept_stat__" "stepsize__"
[15] "treedepth__" "n_leapfrog__"
[17] "divergent__" "energy__" quinn.draw <- quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE)
quinn.draw# A tibble: 10,500 × 5
# Groups: .variable [7]
.chain .iteration .draw .variable .value
<int> <int> <int> <chr> <dbl>
1 1 1 1 b_fSEASONSummer 1.37
2 1 2 2 b_fSEASONSummer 1.14
3 1 3 3 b_fSEASONSummer 1.01
4 1 4 4 b_fSEASONSummer 1.73
5 1 5 5 b_fSEASONSummer 1.57
6 1 6 6 b_fSEASONSummer 1.56
7 1 7 7 b_fSEASONSummer 1.17
8 1 8 8 b_fSEASONSummer 1.26
9 1 9 9 b_fSEASONSummer 0.785
10 1 10 10 b_fSEASONSummer 1.28
# ℹ 10,490 more rowsquinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
summarise_draws(
median,
~ HDInterval::hdi(.x, credMass = 0.95),
rhat,
length,
ess_bulk, ess_tail
)# A tibble: 7 × 9
# Groups: .variable [7]
.variable variable median lower upper rhat length ess_bulk ess_tail
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_DENSITYLow .value 0.950 0.478 1.46 1.00 1500 1502. 1440.
2 b_fSEASONAutumn .value 1.71 0.952 2.64 1.00 1500 1487. 1501.
3 b_fSEASONAutumn:DE… .value 1.02 0.446 2.02 1.00 1500 1473. 1542.
4 b_fSEASONSummer .value 4.13 2.28 6.29 1.00 1500 1410. 1461.
5 b_fSEASONSummer:DE… .value 0.507 0.229 0.957 1.000 1500 1408. 1439.
6 b_fSEASONWinter .value 0.514 0.231 0.834 1.00 1500 1510. 1499.
7 b_fSEASONWinter:DE… .value 0.546 0.168 1.21 0.999 1500 1552. 1385.exceedP <- function(x, Val = 0) mean(x > Val)
quinn.brmsNB |>
tidy_draws() |>
exp() |>
dplyr::select(starts_with("b_")) |>
summarise_draws(
median,
~ HDInterval::hdi(.x, credMass = 0.9),
rhat,
ess_bulk, ess_tail,
~ exceedP(.x, 1)
)# A tibble: 8 × 8
variable median lower upper rhat ess_bulk ess_tail `~exceedP(.x, 1)`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 11.3 8.07 14.9 1.000 1446. 1336. 1
2 b_fSEASONSummer 4.13 2.48 5.80 1.00 1380. 1452. 1
3 b_fSEASONAutumn 1.71 1.05 2.43 0.999 1469. 1489. 0.983
4 b_fSEASONWinter 0.514 0.291 0.783 1.000 1503. 1494. 0.014
5 b_DENSITYLow 0.950 0.585 1.39 1.000 1455. 1393. 0.412
6 b_fSEASONSummer… 0.507 0.262 0.856 1.00 1285. 1442. 0.0347
7 b_fSEASONAutumn… 1.02 0.446 1.73 1.00 1463. 1508. 0.523
8 b_fSEASONWinter… 0.546 0.169 1.02 1.000 1546. 1354. 0.105 We can then summarise this
# A tibble: 7 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_DENSITYLow -0.0513 -0.628 0.424 0.95 median hdci
2 b_fSEASONAutumn 0.537 0.0446 1.01 0.95 median hdci
3 b_fSEASONAutumn:DENSITYLow 0.0206 -0.699 0.744 0.95 median hdci
4 b_fSEASONSummer 1.42 0.931 1.91 0.95 median hdci
5 b_fSEASONSummer:DENSITYLow -0.679 -1.34 0.0263 0.95 median hdci
6 b_fSEASONWinter -0.666 -1.22 -0.0439 0.95 median hdci
7 b_fSEASONWinter:DENSITYLow -0.606 -1.52 0.319 0.95 median hdci We could alternatively express the parameters on the response scale.
# A tibble: 7 × 7
.variable `exp(.value)` .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_DENSITYLow 0.950 0.523 1.51 0.95 median hdci
2 b_fSEASONAutumn 1.71 0.952 2.64 0.95 median hdci
3 b_fSEASONAutumn:DENSITYLow 1.02 0.435 2.01 0.95 median hdci
4 b_fSEASONSummer 4.13 2.31 6.33 0.95 median hdci
5 b_fSEASONSummer:DENSITYLow 0.507 0.229 0.957 0.95 median hdci
6 b_fSEASONWinter 0.514 0.240 0.846 0.95 median hdci
7 b_fSEASONWinter:DENSITYLow 0.546 0.179 1.22 0.95 median hdci quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()
We could alternatively express the parameters on the response scale.
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
group_by(.variable) |>
mutate(.value = exp(.value)) |>
median_hdci()# A tibble: 7 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_DENSITYLow 0.950 0.523 1.51 0.95 median hdci
2 b_fSEASONAutumn 1.71 0.952 2.64 0.95 median hdci
3 b_fSEASONAutumn:DENSITYLow 1.02 0.435 2.01 0.95 median hdci
4 b_fSEASONSummer 4.13 2.31 6.33 0.95 median hdci
5 b_fSEASONSummer:DENSITYLow 0.507 0.229 0.957 0.95 median hdci
6 b_fSEASONWinter 0.514 0.240 0.846 0.95 median hdci
7 b_fSEASONWinter:DENSITYLow 0.546 0.179 1.22 0.95 median hdci quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
ggplot() +
geom_vline(xintercept = 1, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_x_continuous("", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
theme_classic()
Conclusions:
- the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is -0.05. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 0.95.
- the estimated effect of ALG2 vs ALG1 is 0.54 (median) with a standard error of 0.04. The 95% credibility intervals indicate that we are 95% confident that the effect is between 1.01 and 0.95 - e.g. there is a significant positive effect. When back-transformed onto the response scale, we see that barnacle recruitment on ALG2 is 1.71 times higher than that on ALG1. This represents a 71% increase in barnacle recruitment.
- the estimated effect of NB and S are 0.02 and 1.42 respectively, which equate to 0.98 and 0.24 fold declines respectively.
- Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.
## Link scale
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)
## Fractional scale
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
mutate(.value = exp(.value)) |>
ggplot() +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
scale_x_continuous(trans = scales::log2_trans())
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
geom_vline(xintercept = 0, linetype = "dashed")
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = exp(.value), y = .variable)) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans())
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
geom_vline(xintercept = 0, linetype = "dashed")Picking joint bandwidth of 0.0659
## Or on a fractional scale
quinn.brmsNB |>
gather_draws(`.Intercept.*|.*fSEASON.*|.*DENSITY.*`, regex = TRUE) |>
ggplot() +
geom_density_ridges_gradient(
aes(
x = exp(.value),
y = .variable,
fill = stat(x)
),
alpha = 0.4, colour = "white",
quantile_lines = TRUE,
quantiles = c(0.025, 0.975)
) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans()) +
scale_fill_viridis_c(option = "C")Picking joint bandwidth of 0.0951
This is purely a graphical depiction on the posteriors.
# A tibble: 1,500 × 21
.chain .iteration .draw b_Intercept b_fSEASONSummer b_fSEASONAutumn
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 2.30 1.37 0.979
2 1 2 2 2.77 1.14 0.0894
3 1 3 3 2.79 1.01 0.380
4 1 4 4 2.14 1.73 0.922
5 1 5 5 2.21 1.57 0.697
6 1 6 6 2.21 1.56 0.963
7 1 7 7 2.51 1.17 0.572
8 1 8 8 2.75 1.26 0.625
9 1 9 9 2.82 0.785 0.754
10 1 10 10 2.53 1.28 0.509
# ℹ 1,490 more rows
# ℹ 15 more variables: b_fSEASONWinter <dbl>, b_DENSITYLow <dbl>,
# `b_fSEASONSummer:DENSITYLow` <dbl>, `b_fSEASONAutumn:DENSITYLow` <dbl>,
# `b_fSEASONWinter:DENSITYLow` <dbl>, shape <dbl>, Intercept <dbl>,
# lprior <dbl>, lp__ <dbl>, accept_stat__ <dbl>, stepsize__ <dbl>,
# treedepth__ <dbl>, n_leapfrog__ <dbl>, divergent__ <dbl>, energy__ <dbl># A tibble: 1,500 × 10
.chain .iteration .draw b_fSEASONSummer b_fSEASONAutumn b_fSEASONWinter
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 1.37 0.979 -0.362
2 1 2 2 1.14 0.0894 -0.959
3 1 3 3 1.01 0.380 -0.558
4 1 4 4 1.73 0.922 -0.469
5 1 5 5 1.57 0.697 -0.269
6 1 6 6 1.56 0.963 -0.682
7 1 7 7 1.17 0.572 -0.742
8 1 8 8 1.26 0.625 -0.723
9 1 9 9 0.785 0.754 -1.13
10 1 10 10 1.28 0.509 -0.708
# ℹ 1,490 more rows
# ℹ 4 more variables: b_DENSITYLow <dbl>, `b_fSEASONSummer:DENSITYLow` <dbl>,
# `b_fSEASONAutumn:DENSITYLow` <dbl>, `b_fSEASONWinter:DENSITYLow` <dbl>Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.# A tibble: 1,500 × 12
b_Intercept b_fSEASONSummer b_fSEASONAutumn b_fSEASONWinter b_DENSITYLow
<dbl> <dbl> <dbl> <dbl> <dbl>
1 2.30 1.37 0.979 -0.362 -0.111
2 2.77 1.14 0.0894 -0.959 -0.703
3 2.79 1.01 0.380 -0.558 -0.626
4 2.14 1.73 0.922 -0.469 0.414
5 2.21 1.57 0.697 -0.269 0.198
6 2.21 1.56 0.963 -0.682 0.266
7 2.51 1.17 0.572 -0.742 -0.209
8 2.75 1.26 0.625 -0.723 -0.324
9 2.82 0.785 0.754 -1.13 -0.596
10 2.53 1.28 0.509 -0.708 0.0113
# ℹ 1,490 more rows
# ℹ 7 more variables: `b_fSEASONSummer:DENSITYLow` <dbl>,
# `b_fSEASONAutumn:DENSITYLow` <dbl>, `b_fSEASONWinter:DENSITYLow` <dbl>,
# shape <dbl>, Intercept <dbl>, lprior <dbl>, lp__ <dbl>Region of Practical Equivalence
[1] 0.2754809Possible multicollinearity between b_fSEASONSummer:DENSITYLow and
b_DENSITYLow (r = 0.71). This might lead to inappropriate results. See
'Details' in '?rope'.# Proportion of samples inside the ROPE [-0.28, 0.28]:
Parameter | inside ROPE
--------------------------------------
Intercept | 0.00 %
fSEASONSummer | 0.00 %
fSEASONAutumn | 11.94 %
fSEASONWinter | 6.81 %
DENSITYLow | 77.18 %
fSEASONSummer:DENSITYLow | 12.92 %
fSEASONAutumn:DENSITYLow | 58.57 %
fSEASONWinter:DENSITYLow | 22.82 %Possible multicollinearity between b_fSEASONSummer:DENSITYLow and
b_DENSITYLow (r = 0.71). This might lead to inappropriate results. See
'Details' in '?rope'.
## Or based on fractional scale
quinn.brmsNB |>
as_draws_df("^b_fSEASON.*|^b_DENSITY.*", regex = TRUE) |>
exp() |>
## equivalence_test(range = c(0.755, 1.32))
rope(range = c(0.755, 1.32))# Proportion of samples inside the ROPE [0.76, 1.32]:
Parameter | inside ROPE
--------------------------------------
fSEASONSummer | 0.00 %
fSEASONAutumn | 11.73 %
fSEASONWinter | 6.81 %
DENSITYLow | 76.90 %
fSEASONSummer:DENSITYLow | 12.92 %
fSEASONAutumn:DENSITYLow | 58.43 %
fSEASONWinter:DENSITYLow | 22.89 %quinn.mcmc <-
quinn.brmsNB |>
as_draws_df("^b_fSEASON.*|^b_DENSITY.*", regex = TRUE) |>
exp()
quinn.mcmc |>
rope(range = c(0.755, 1.32))# Proportion of samples inside the ROPE [0.76, 1.32]:
Parameter | inside ROPE
--------------------------------------
fSEASONSummer | 0.00 %
fSEASONAutumn | 11.73 %
fSEASONWinter | 6.81 %
DENSITYLow | 76.90 %
fSEASONSummer:DENSITYLow | 12.92 %
fSEASONAutumn:DENSITYLow | 58.43 %
fSEASONWinter:DENSITYLow | 22.89 %## note, the following is not quit correct, it does not get the CI correct
quinn.mcmc |>
rope(range = c(0.755, 1.32)) |>
plot(data = quinn.brmsNB)
# Test for Practical Equivalence
ROPE: [0.76 1.32]
Parameter | H0 | inside ROPE | 95% HDI
-----------------------------------------------------------------
fSEASONSummer | Rejected | 0.00 % | [2.48, 6.59]
fSEASONAutumn | Undecided | 11.73 % | [1.05, 2.77]
fSEASONWinter | Undecided | 6.81 % | [0.28, 0.92]
DENSITYLow | Undecided | 76.90 % | [0.55, 1.59]
fSEASONSummer:DENSITYLow | Undecided | 12.92 % | [0.26, 1.05]
fSEASONAutumn:DENSITYLow | Undecided | 58.43 % | [0.50, 2.12]
fSEASONWinter:DENSITYLow | Undecided | 22.89 % | [0.21, 1.35]quinn.brmsNB |> modelsummary(
statistic = c("conf.low", "conf.high"),
shape = term ~ statistic,
exponentiate = TRUE
)Warning:
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
This warning appears once per session.| (1) | |||
|---|---|---|---|
| Est. | 2.5 % | 97.5 % | |
| b_Intercept | 11.259 | 7.868 | 16.737 |
| b_fSEASONSummer | 4.127 | 2.478 | 6.590 |
| b_fSEASONAutumn | 1.710 | 1.052 | 2.765 |
| b_fSEASONWinter | 0.514 | 0.277 | 0.922 |
| b_DENSITYLow | 0.950 | 0.547 | 1.591 |
| b_fSEASONSummer × DENSITYLow | 0.507 | 0.265 | 1.048 |
| b_fSEASONAutumn × DENSITYLow | 1.021 | 0.500 | 2.119 |
| b_fSEASONWinter × DENSITYLow | 0.546 | 0.215 | 1.354 |
| Num.Obs. | 42 | ||
| R2 | 0.732 | ||
| ELPD | -148.3 | ||
| ELPD s.e. | 5.9 | ||
| LOOIC | 296.7 | ||
| LOOIC s.e. | 11.8 | ||
| WAIC | 295.7 | ||
| RMSE | 7.48 | ||
11 Further investigations
fSEASON = Spring:
contrast ratio lower.HPD upper.HPD
High / Low 0.898 0.341 1.67
fSEASON = Summer:
contrast ratio lower.HPD upper.HPD
High / Low 2.187 0.995 3.90
fSEASON = Autumn:
contrast ratio lower.HPD upper.HPD
High / Low 1.069 0.417 2.01
fSEASON = Winter:
contrast ratio lower.HPD upper.HPD
High / Low 2.122 0.475 5.43
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95 ## absolute response scale
quinn.rstanarmNB |>
emmeans(~ DENSITY | fSEASON, type = "link") |>
regrid() |>
pairs()fSEASON = Spring:
contrast estimate lower.HPD upper.HPD
High - Low -1.14 -9.41 7.22
fSEASON = Summer:
contrast estimate lower.HPD upper.HPD
High - Low 26.30 1.54 53.14
fSEASON = Autumn:
contrast estimate lower.HPD upper.HPD
High - Low 1.28 -14.93 16.27
fSEASON = Winter:
contrast estimate lower.HPD upper.HPD
High - Low 2.91 -1.10 7.68
Point estimate displayed: median
HPD interval probability: 0.95 quinn.em <- quinn.rstanarmNB |>
emmeans(~ DENSITY | fSEASON, type = "link") |>
pairs() |>
gather_emmeans_draws() |>
mutate(Fit = exp(.value))
head(quinn.em)# A tibble: 6 × 7
# Groups: contrast, fSEASON [1]
contrast fSEASON .chain .iteration .draw .value Fit
<fct> <fct> <int> <int> <int> <dbl> <dbl>
1 High - Low Spring NA NA 1 0.0120 1.01
2 High - Low Spring NA NA 2 -0.439 0.645
3 High - Low Spring NA NA 3 -0.154 0.858
4 High - Low Spring NA NA 4 0.0245 1.02
5 High - Low Spring NA NA 5 0.106 1.11
6 High - Low Spring NA NA 6 -0.0619 0.940g2 <- quinn.em |>
group_by(contrast, fSEASON) |>
median_hdci() |>
ggplot() +
geom_vline(xintercept = 1, linetype = "dashed") +
geom_pointrange(aes(x = Fit, y = fSEASON, xmin = Fit.lower, xmax = Fit.upper)) +
scale_x_continuous("Effect size (High/Low)", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
theme_classic()
g2
ggplot(quinn.em, aes(x = Fit)) +
geom_histogram() +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous("Effect size (High/Low)", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
facet_wrap(fSEASON ~ contrast, scales = "free")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# A tibble: 4 × 8
contrast fSEASON Fit .lower .upper .width .point .interval
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 High - Low Spring 0.898 0.394 1.74 0.95 median hdci
2 High - Low Summer 2.19 0.995 3.90 0.95 median hdci
3 High - Low Autumn 1.07 0.440 2.03 0.95 median hdci
4 High - Low Winter 2.12 0.475 5.43 0.95 median hdci `summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.# A tibble: 4 × 3
# Groups: contrast [1]
contrast fSEASON P
<fct> <fct> <dbl>
1 High - Low Spring 0.388
2 High - Low Summer 0.992
3 High - Low Autumn 0.571
4 High - Low Winter 0.919## Probability of effect greater than 10%
quinn.em |>
group_by(contrast, fSEASON) |>
summarize(P = mean(Fit > 1.1))`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.# A tibble: 4 × 3
# Groups: contrast [1]
contrast fSEASON P
<fct> <fct> <dbl>
1 High - Low Spring 0.287
2 High - Low Summer 0.977
3 High - Low Autumn 0.467
4 High - Low Winter 0.891## Using summarise_draws
quinn.rstanarmNB |>
emmeans(~ DENSITY | fSEASON, type = "link") |>
pairs() |>
tidy_draws() |>
exp() |>
summarise_draws(median,
HDInterval::hdi,
P = ~ mean(.x > 1)
)# A tibble: 4 × 5
variable median lower upper P
<chr> <dbl> <dbl> <dbl> <dbl>
1 contrast High - Low fSEASON Spring 0.898 0.341 1.67 0.388
2 contrast High - Low fSEASON Summer 2.19 0.995 3.90 0.992
3 contrast High - Low fSEASON Autumn 1.07 0.417 2.01 0.571
4 contrast High - Low fSEASON Winter 2.12 0.475 5.43 0.919newdata <- with(quinn, expand.grid(
fSEASON = levels(fSEASON),
DENSITY = levels(DENSITY)
))
Xmat <- model.matrix(~ fSEASON * DENSITY, data = newdata)
as.matrix(quinn.rstanarmNB) |> head() parameters
iterations (Intercept) fSEASONSummer fSEASONAutumn fSEASONWinter DENSITYLow
[1,] 2.331955 1.455598 0.7542465 -0.4274555 -0.01198242
[2,] 2.257843 1.392760 0.8617720 -0.3243542 0.43857257
[3,] 1.944257 2.055757 1.0911906 -0.3493536 0.15370757
[4,] 2.399954 1.274637 0.6597981 -0.8023593 -0.02454189
[5,] 2.237445 1.526110 1.2116360 -0.7028590 -0.10646826
[6,] 2.261734 1.585028 0.6853106 -0.4497903 0.06193139
parameters
iterations fSEASONSummer:DENSITYLow fSEASONAutumn:DENSITYLow
[1,] -0.8273987 0.02463655
[2,] -1.2303391 -0.65127644
[3,] -0.8041500 -0.31550231
[4,] -0.4068205 0.16244078
[5,] -0.4684958 -0.48006959
[6,] -0.7467988 0.07305836
parameters
iterations fSEASONWinter:DENSITYLow reciprocal_dispersion
[1,] -0.3375743 4.427056
[2,] -2.0220195 4.663400
[3,] -1.3351653 4.475706
[4,] -0.3388785 2.882597
[5,] -0.4987464 3.866986
[6,] -0.9840772 4.305455## coefs <- as.matrix(quinn.rstanarmNB)
coefs <- as.matrix(as.data.frame(quinn.rstanarmNB) |>
dplyr:::select(-reciprocal_dispersion)) |>
as.matrix()
fit <- exp(coefs %*% t(Xmat))
newdata <- newdata |>
cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
head(newdata) fSEASON DENSITY term estimate std.error conf.low conf.high
1 Spring High 1 10.431107 2.728881 5.816823 16.155118
2 Summer High 2 50.693111 12.082408 28.431714 73.857594
3 Autumn High 3 20.646629 5.138216 12.640919 31.889828
4 Winter High 4 5.886958 1.708075 3.051383 9.209016
5 Spring Low 5 11.723876 3.280842 6.174902 18.344092
6 Summer Low 6 23.126044 5.735582 13.086952 34.588108ggplot(newdata, aes(y = estimate, x = fSEASON, fill = DENSITY)) +
geom_blank() +
geom_line(aes(x = as.numeric(fSEASON), ymin = conf.low, ymax = conf.high, linetype = DENSITY),
position = position_dodge(0.2)
) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
shape = 21,
position = position_dodge(0.2)
)Warning in geom_line(aes(x = as.numeric(fSEASON), ymin = conf.low, ymax =
conf.high, : Ignoring unknown aesthetics: ymin and ymax
# Compare high and low in each season
# via contrasts
newdata <- with(quinn, expand.grid(
fSEASON = levels(fSEASON),
DENSITY = levels(DENSITY)
))
## factor differences
Xmat <- model.matrix(~ fSEASON * DENSITY, data = newdata)
Xmat.high <- Xmat[newdata$DENSITY == "High", ]
Xmat.low <- Xmat[newdata$DENSITY == "Low", ]
Xmat.density <- Xmat.high - Xmat.low
rownames(Xmat.density) <- levels(quinn$fSEASON)
coefs <- as.matrix(as.data.frame(quinn.rstanarmNB) |> dplyr:::select(-reciprocal_dispersion))
fit <- exp(coefs %*% t(Xmat.density))
tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")# A tibble: 4 × 5
term estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 Spring 0.956 0.365 0.341 1.67
2 Summer 2.32 0.790 0.995 3.90
3 Autumn 1.14 0.450 0.417 2.01
4 Winter 2.52 1.56 0.475 5.43## or absolute
fit.high <- coefs %*% t(Xmat.high)
fit.low <- coefs %*% t(Xmat.low)
fit <- exp(fit.high) - exp(fit.low)
# fit = exp(fit.high - fit.low)
tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")# A tibble: 4 × 5
term estimate std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 1 -1.29 4.24 -9.41 7.22
2 2 27.6 13.4 1.54 53.1
3 3 0.977 7.94 -14.9 16.3
4 4 2.95 2.24 -1.10 7.68Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.fSEASON = Spring:
contrast ratio lower.HPD upper.HPD
High / Low 1.05 0.555 1.68
fSEASON = Summer:
contrast ratio lower.HPD upper.HPD
High / Low 2.04 1.198 3.22
fSEASON = Autumn:
contrast ratio lower.HPD upper.HPD
High / Low 1.04 0.510 1.75
fSEASON = Winter:
contrast ratio lower.HPD upper.HPD
High / Low 1.94 0.671 4.13
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95 ## absolute response scale
quinn.brmsNB |>
emmeans(~ DENSITY | fSEASON, type = "link") |>
regrid() |>
pairs()Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.fSEASON = Spring:
contrast estimate lower.HPD upper.HPD
High - Low 0.566 -5.223 6.74
fSEASON = Summer:
contrast estimate lower.HPD upper.HPD
High - Low 23.671 6.579 41.87
fSEASON = Autumn:
contrast estimate lower.HPD upper.HPD
High - Low 0.727 -12.333 12.30
fSEASON = Winter:
contrast estimate lower.HPD upper.HPD
High - Low 2.726 -0.884 6.30
Point estimate displayed: median
HPD interval probability: 0.95 quinn.em <- quinn.brmsNB |>
emmeans(~ DENSITY | fSEASON, type = "link") |>
pairs() |>
gather_emmeans_draws() |>
mutate(Fit = exp(.value))Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 6 × 7
# Groups: contrast, fSEASON [1]
contrast fSEASON .chain .iteration .draw .value Fit
<fct> <fct> <int> <int> <int> <dbl> <dbl>
1 High - Low Spring NA NA 1 0.111 1.12
2 High - Low Spring NA NA 2 0.703 2.02
3 High - Low Spring NA NA 3 0.626 1.87
4 High - Low Spring NA NA 4 -0.414 0.661
5 High - Low Spring NA NA 5 -0.198 0.821
6 High - Low Spring NA NA 6 -0.266 0.766g2 <- quinn.em |>
group_by(contrast, fSEASON) |>
median_hdci() |>
ggplot() +
geom_vline(xintercept = 1, linetype = "dashed") +
geom_pointrange(aes(x = Fit, y = fSEASON, xmin = Fit.lower, xmax = Fit.upper)) +
scale_x_continuous("Effect size (High/Low)", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
theme_classic()
g2
ggplot(quinn.em, aes(x = Fit)) +
geom_histogram() +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous("Effect size (High/Low)", trans = scales::log2_trans(), breaks = unique(as.vector(2^(0:4 %o% c(-1, 1))))) +
facet_wrap(fSEASON ~ contrast, scales = "free")`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# A tibble: 4 × 8
contrast fSEASON Fit .lower .upper .width .point .interval
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 High - Low Spring 1.05 0.578 1.71 0.95 median hdci
2 High - Low Summer 2.04 1.20 3.22 0.95 median hdci
3 High - Low Autumn 1.04 0.510 1.75 0.95 median hdci
4 High - Low Winter 1.94 0.671 4.13 0.95 median hdci `summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.# A tibble: 4 × 3
# Groups: contrast [1]
contrast fSEASON P
<fct> <fct> <dbl>
1 High - Low Spring 0.588
2 High - Low Summer 0.997
3 High - Low Autumn 0.553
4 High - Low Winter 0.935## Probability of effect greater than 10%
quinn.em |>
group_by(contrast, fSEASON) |>
summarize(P = mean(Fit > 1.1))`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.# A tibble: 4 × 3
# Groups: contrast [1]
contrast fSEASON P
<fct> <fct> <dbl>
1 High - Low Spring 0.425
2 High - Low Summer 0.989
3 High - Low Autumn 0.423
4 High - Low Winter 0.90312 Summary figures
newdata <- quinn.rstanarmNB |>
emmeans(~ fSEASON | DENSITY, type = "response") |>
as.data.frame()
head(newdata) fSEASON DENSITY prob lower.HPD upper.HPD
Spring High 10.08644 5.816823 16.15512
Summer High 49.32256 28.431714 73.85759
Autumn High 19.85178 12.640919 31.88983
Winter High 5.65834 3.051383 9.20902
Spring Low 11.29123 6.174902 18.34409
Summer Low 22.42250 13.086952 34.58811
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95 g1 <- ggplot(newdata, aes(y = prob, x = fSEASON, color = DENSITY)) +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD),
position = position_dodge(width = 0.2)
) +
theme_classic()
g1 + g2
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead. fSEASON DENSITY prob lower.HPD upper.HPD
Spring High 11.25880 7.629901 16.05913
Summer High 46.44690 31.426065 63.52073
Autumn High 19.34153 13.173481 27.28439
Winter High 5.72578 3.417179 8.89959
Spring Low 10.63526 6.703797 15.44913
Summer Low 22.46584 15.270754 32.59736
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95 g1 <- ggplot(newdata, aes(y = prob, x = fSEASON, color = DENSITY)) +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD),
position = position_dodge(width = 0.2)
) +
theme_classic()
g1 + g2
13 Observation-level random effects
13.1 brms
quinn <- quinn |>
group_by(fSEASON, DENSITY) |>
mutate(Obs = factor(1:n()))
quinn.form <- bf(RECRUITS ~ fSEASON * DENSITY + (1 | Obs), family = poisson(link = "log"))
get_prior(quinn.form, data = quinn) prior class coef group resp dpar
(flat) b
(flat) b DENSITYLow
(flat) b fSEASONAutumn
(flat) b fSEASONAutumn:DENSITYLow
(flat) b fSEASONSummer
(flat) b fSEASONSummer:DENSITYLow
(flat) b fSEASONWinter
(flat) b fSEASONWinter:DENSITYLow
student_t(3, 2.6, 2.5) Intercept
student_t(3, 0, 2.5) sd
student_t(3, 0, 2.5) sd Obs
student_t(3, 0, 2.5) sd Intercept Obs
nlpar lb ub source
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
0 default
0 (vectorized)
0 (vectorized)quinn.brmsU <- brm(quinn.form,
data = quinn,
refresh = 0,
chains = 3,
iter = 5000,
thin = 5,
warmup = 2000
)Compiling Stan program...Start samplingWarning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problemsWarning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.quinn.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = quinn$RECRUITS,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(quinn.resids)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.DENSITY = High:
fSEASON rate lower.HPD upper.HPD
Spring 9.72548 6.28423 13.38155
Summer 46.88011 32.28727 60.64975
Autumn 19.29252 12.88027 25.56589
Winter 5.50748 3.55364 8.26443
DENSITY = Low:
fSEASON rate lower.HPD upper.HPD
Spring 10.69229 6.95720 14.79015
Summer 21.58560 15.50071 29.06059
Autumn 16.39288 10.90293 22.47985
Winter 2.38235 0.87857 4.49520
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95 ggplot(newdata, aes(y = rate, x = fSEASON, color = DENSITY)) +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD),
position = position_dodge(width = 0.2)
)
14 References
Quinn, G. P. 1988. “Ecology of the Intertidal Pulmonate Limpet Siphonaria Diemenensis Quoy Et Gaimard. II Reproductive Patterns and Energetics.” Journalofexperimentalmarinebiologyandecology 117: 137–56.